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I.  INTRODUCTION 


This  report  describes  the  use  of  computer  simulation  tech- 
niques to  carry  out  ntunerical  modelling  of  dynamical  properties 
associated  with  ionospheric  and  plasmaspheric  boundaries.  Two 
specific  cases  are  dealt  with:  (1)  changes  in  the  equatorial 
plasmapause  location  during  periods  of  geomagnetic  activity,  and 
C2)  the  effects  of  vertical  plasma  diffusion  upon  aurorally- 
induced  ionization  enhancements  at  ionospheric  heights.  For  the 
former  case,  a scheme  was  developed  whereby  changes  in  magneto- 
spheric  convection  patterns  could  be  related  to  plasmapause 
variations  by  using  computer  generated  "maps"  of  an  ensemble  of 
equatorial  "test  particles"  to  define  equilibrium  and  non-equi- 
librium configurations  of  the  plasmapause.  For  the  latter  case, 
finite  element  simulation  (FES)  techniques  were  used  to  construct 
a one- dimensional  model  of  the  temporal  evolution  of  ionospheric 
plasma  produced  by  energetic  particles  as  found,  for  example, 
beneath  the  magnetospheric  cusps  or  near  the  poleward  wall  of 
the  F-region  trough.  In  both  cases,  the  results  obtained  show 
that  plasma  transport  effects  either  dominate  or  at  least  signif- 
icantly contribute  to  the  spatial  and/or  temporal  characteristics 
of  the  feature  in  question. 

In  Chapter  2,  our  magnetospheric-convection-dominated  model 
for  the  formation  of  the  plasmasphere  is  described.  Calculations 
of  convection  effects  in  the  equatorial  plane  are  used  to  address 
two  aspects  of  time- dependent  plasmapause  dynamics:  (1)  the 
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Transition  Time  from  one  stable,  equilibrium  plasmapause  config- 
uration to  another,  and  (2)  Disturbance  Scenarios,  i.e.,  simu- 
lations of  the  hour-by-hour  change  in  the  plasmapause  location 
on  the  nightside  associated  with  specific  variations  in  geomag- 
netic activity. 

In  Chapter  3,  we  describe  the  development  of  a computational 
scheme  based  on  the  Finite  Element  Simulation  technique  for 
handling  particle-produced  ionization  at  high  latitudes.  Initial 
results  obtained  by  applying  a one-dimensional  version  of  this 
computational  scheme  to  the  cusp-related  electron  density  en- 
hancements on  the  dayside  are  also  presented. 
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2.  SIMULATION  OF  PLASMAPAUSE  DYNAMICS 

2.1  Background  and  Physical  Assumptions 

Most  of  the  fundamental  theoretical  studies  of  the  origin 
and  behavior  of  the  equatorial  plasmapause  (e.g.,  Nishida,  1966; 
Brice,  1967;  Kavanagh  et  al. , 1968)  deal  with  the  application  of 
magnetospheric  convection  patterns  originally  suggested  by  Axford 
and  Hines  C1961)  to  the  classic  plasmapause  geometry  reported  by 
Carpenter  tl966) . The  notion  of  a contracting  plasmapause  due 
to  enhanced  convection  patterns  during  periods  of  increased  geo- 
magnetic activity  is  now  a widely  accepted  view,  and  is,  in  fact, 
the  scheme  generally  said  to  explain  why  statistical  results  for 
the  plasmapause  location  (L__)  vs.  an  index  like  Kp  seem  to  pro- 
vide  reasonable  descriptions  of  average  behavior  (Ry croft  and 
Thomas,  1970;  Kohnlein  and  Raitt,  1977).  Yet,  these  simple 
statistical  patterns  of  L vs.  Kp  are  not  expected  to  provide 

Crir 

an  actual  description  of  the  time  dependence  of  plasmapause 
changes  during  a specific  storm's  Kp-scenario.  The  real  plasma- 
pause cannot  adjust  instantaneously  to  a new  configuration  for 
each  Kp  step,  and  within  a given  Kp  interval,  the  L position 

Pc' 

must  be  a function  of  previous  magnetic  history. 

The  concept  of  continuously  monitoring  the  plasmapause  loca- 
tion via  ground-based  or  satellite  techniques  is  not  a very 
realistic  possibility.  Whistler  measurements  from  a fixed  site 
sample  a family  of  field  lines,  but  the  measurements  are  obvi- 
ously fixed  to  the  Earth's  rotating  frame  of  reference  and  there- 
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fore  only  one  local  time  sector  is  observed  at  a time.  Satellite 
measurements  in  or  near  the  equatorial  plane  offer  the  best 
method  of  determining  Lpp  as  a function  of  local  time;  due  to  the 
need  to  sample  both  in  the  radial  and  azimuthal  directions,  a 
single  equatorial  satellite  does  not  readily  provide  good  L 

PP 

time  histories.  By  far  the  easiest  scheme  would  be  to  determine 

the  plasmapause  signature  at  ionospheric  heights  where  several 

high  inclination  satellites  could  provide  L vs  local  time 

PP 

morphologies.  However,  specifying  the  plasmapause  signature  in 

the  F-region  is  not  a widely  agreed  upon  concept.  At  various 

times,  clear  L signatiares  in  the  topside  ionosphere  appear  as 

+ + 

(1)  an  abrupt  increase  in  the  H /O  transition  height  (i.e.,  the 
"light  ion  trough'")  , (2)  a plasma  temperature  peak  or  (3)  a steep 
total  ion  gradient  (i.e.,  an  electron  density  ledge).  The  vari- 
ous altitude  and  latitude  characteristics  of  the  quiet-time 
trough  have  recently  been  reviewed  by  Mendillo  and  Chacko  (1977) . 

In  testing  a numerical  simulation  scheme  for  specific  storm-time 
scenarios  (see  section  2.2),  we  decided  to  rely  upon  characteristic 
C3)  above,  that  is,  we  attempted  to  reproduce  satellite  probe 
determinations  of  the  trough/plasmapause  location  in  the  topside 
F-region  (h  = 800  - 3200  km) . The  studies  carried  out  by 
Bewersdorff  and  Sagalyn  (1972)  and  Wildman  et  al.  (1976)  support 
the  view  that  d\iring  periods  of  increased  geomagnetic  activity, 
the  steep  equatorward  edge  of  the  topside  plasma  trough  is  prob- 
ably a good  indicator  of  the  plasmapause  location. 
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In  developing  a model  for  plasmaspheric  dynamics,  we  assume 
a completely  convection-dominated  scheme.  By  neglecting  iono- 
spheric replenishment  of  the  plasmasphere , we  limit  our  results 
to  nighttime  effects  and  to  simulation  times  not  much  longer  than 
the  co-rotation  period  (24  hours) . Both  constraints  again  seem 
reasonable  for  periods  of  increasing  geomagnetic  activity. 

The  magnetospheric  convection  code  developed  (called  Program 
MAGCON)  relies  upon  computer-generated  (i.e.,  CALCOMP)  "maps"  of 
an  ensemble  of  plasmaspheric  "test  particles"  to  define  equilib- 
rium and  transitional  configurations  of  the  plasmapause.  The 
model  was  used  to  address  two  specific  aspects  of  plasmaspheric 
dynamics : 

(1)  The  Transition  Time  (Tp)  from  one  stable,  equilibrium 
plasmapause  configuration  to  another. 

C2)  Disturbance  Scenario^,  i.e.,  the  hour-by-hour  changes 

in  the  plasmapause  location  at  a given  local  time  (I'pp(LT)) 
associated  with  specified  variations  in  geomagnetic  activity. 
All  of  the  model  studies  carried  out  used  the  geometry  and 
electric  field  configurations  depicted  in  Figures  1 and  2 
CMendillo,  1973;  Mendillo  and  Papagiannis,  1971).  For  any  given 
test  particle  at  a specified  local  time  and  L-value,  the  instan- 
taneous motion  of  the  particle  is  assumed  to  be  given  by 
2 

W = (E  X b)/B  , where  E is  the  resultant  electric  field  due  to 
co-rotation  directed  radially  inward)  and  convection 

directed  from  dawn  to  dusk) , and  B is  the  dipole  field  strength 
in  the  equatorial  plane.  A "fluid"  of  particles  is  followed  by 
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Figtire  1.  Equatorial  Cross-Section  of  the  Magnetosphere 
(9  * 90”).  Concentric  ri^s  unit  L-shells,  arrows  in 
the  -X  direction  give  the  E x B flow  lines  resulting  from 
the  dawn-dusk  electric  field  y 
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Fig.  3.  The  solid  curve  gives  the  radial  be- 
havior of  the  pl.asaapause  at  dusk  (£„)  v»rsi!s 
Kp  from  equations  IS  an'i  .'0  of  the  model  pre- 
sented. The  dashed  curve  represcr's  the  linear 
Ln  dependence  .‘usaested  by  Vasyiiunrs  (lOfiSl. 
The  crosses  show  the  beh.avior  of  the  ionospheric 
trough  at  ISOO  LT  from  Thovnn,^  a-,d  dnaVerrs 
[1063];  the  dors  .are  H*  pbsmapauso  positions 
near  2400  LT  from  Taylor  el  al.  [19691. 
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keeping  track  of  all  the  individual  motions  using  a 10  minute 
time  step.  The  initial  distribution  of  particles  was  taken  to 
be  circular,  with  a particle  at  each  % L-value  from  L = 1.5  to 
5.5,  and  spaced  every  15  minutes  of  local  time.  A total  of  864 
individual  points  are  therefore  available  to  define  subsequent 
transformations . 

2.2.  Simulation  of  Plasmapause  Transition  Times  (Tp) 

For  a given  dawn-to-dusk  electric  field  (assumed  to  be  tem- 
porally and  spatially  constant) , a unique  plasmapause  exists 
CKavanagh  et  al.,  1968;  Chen,  1970)  which  defines  the  boundary 
between  plasma  particles  which  share  in  the  co- rotational  field 
associated  with  the  earth  and  those  which  are  lost  due  to  mag- 
netospheric  convection.  The  shape  of  the  resultant  plasmapause 
for  such  a case  is  the  familiar  "tear-drop"  described  by  Chen 
C1970) . Thus,  our  initial  circular  distribution  of  particles 
becomes  distorted  into  a "teeur-drop"  configuration  due  to  the 

field.  It  typically  tcdces  cdjout  30  hours  of  simulation  time 
for  a stable  configuration  to  develop,  i.e.,  the  remaining  fluid 
particles  (perhaps  50  to  75%  of  the  original  number)  are  trapped 
in  motion  about  the  Earth  in  a region  defined  by  the  appropriate 
theoretical  "tear-drop"  plasmapause  (see  Figure  6 in  next  section) . 
Given  diis  equilibrium  plasmapause,  we  are  interested  in  subse- 
quent transformations  caused  by  a stronger  Ejjjj  field.  Clearly, 
the  new  plasmapause  position  is  known  in  advance  since  it  depends 
only  on  the  new  value  of  E^j^.  The  actual  time  it  takes  to  go 
from  one  "tear-drop"  to  another  is  not  known,  and  thus  program 
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MAGCON  lends  itself  nicely  to  this  type  of  "dynamic  calculation". 

In  an  earlier  work  (Mendillo  and  PapagianniS/  1971)  we 
derived  an  expression  for  the  dependence  of  the  magnetospheric 
electric  field  on  the  solar  wind  speed  and  related  the 

results  to  plasmapause  changes  at  1800  LT  through  the  Kp  index 
Csee  Figure  3).  Using  this  model  for  vs.  Kp,  we  computed  the 
set  of  10  equilibrium  "tear-drop"  plasmapause  configurations  for 
Kp  = 0 to  9 . These  are  plotted  on  LT/L-value  grid  in  Figure  4, 
and  tabulated  values  of  selected  parameters  are  given  in  Table  1. 

As  can  be  seen  from  Figure  4,  the  largest  changes  in  L as 
Kp  increases  occur  in  the  dusk  sector.  While  the  total  rauige  of 
L values  experienced  near  1800  LT  seems  consistent  with  avail- 
able  data ’Csee  Figure  3),  the  rather  low  values' at  00  LT  and 

pp 

06  LT  suggest  that  the  simple  "tear-drop"  model,  while  qualita- 
tively giving  the  correct  dawn  - dusk  asymetry,  does  not  hold  up 
to  quantitative  comparisons  with  real  data  at  all  local  times. 

This  "improper  shape"  of  the  plasmapause  will  be  discussed  in  a 
later  report.  For  the  present  purpose  of  obtaining  estimates  of 
Transition  Times  (Tp)  for  the  plasmapause  configuration  associated 
with  one  Kp  to  that  associated  with  another,  we  decided  to  use 
the  existing  tear-drop  configurations. 

In  operationally  computing  Tp  (say  from  Kp  = 1 to  3) , the 
ensemble  of  particles  occupying  the  Kp  = 1 tear-drop  experiences 
some  "peeling-off"  in  the  12-18  LT  sector  and  compression  effects 
on  the  nightside.  Particles  are  considered  to  be  lost  if  they 


10 


TABLE  1 

EQUILIBRIUM  PLASMAPAUSE  MODEL  NORMALIZED  TO  DUSK  SECTOR 


Kp 

V(sw) 

E 

conv 

^PP 

L 

PP 

L 

PP 

(km/sec) 

(mv/m) 

(kv/re) 

(dusk) 

(noon/midn) 

(dawn) 

0 

330.0 

0.22 

1.39 

8.16 

4.08 

3.  38 

1 

397.5 

0.32 

2.02 

6.77 

3.38 

2.80 

2 

465.0 

0.43 

2.76 

5.79 

2.89 

2.40 

3 

532,6 

0.57 

3.62 

5.05 

2.53 

2.09 

4 

600.1 

0.72 

4.60 

4.48 

2.24 

1.  86 

5 

667.6 

0.89 

5.69 

4.03 

2.02 

1.67 

6 

735.1 

1.08 

6.91 

3.66 

1.  83 

1.52 

7 

802.6 

1.29 

• 8.23 

3.35 

1.68 

1.  39 

8 

870.2 

1.52 

9.67 

3.09 

1.55 

1.28 

9 

937.7 

1.76 

11.23 

2.87 

1.  44 

1.19 

11 


convect  to  an  L-value  greater  than  the  new  L (18  LT)  value,  or 

ir  Jr 

move  to  L = 1.  We  define  "100%  equilibrium"  to  occur  when  liter- 
ally all  the  particles  capable  of  being  lost  are  lost.  Similarly, 
the  25%,  50%,  and  90%  equilibrium  times  are  the  time  intervals  it 
takes  for  25%,  50%,  and  90%  of  the  particles  going  to  be  lost  are 
lost,  respectively.  A set  of  such  Transition  Time  results  are 
displayed  in  Table  2.  Note  that  the  format  gives  the  Tp(25%), 
TpC50%),  TpC90%),  and  Tp(100%)  times  for  transformations  of  Kp 
Cinitial)  to  Kp  (final).  In  examining  the  Tp(l00%)  results,  one 
ceui  see  that  Tp  is  large  (-30  hours)  , independent  of  whether  the 
Kp  step  is  small  (i.e.,  Kp  = 0-*-l,  1-^2,  0->2,  6->-8)  or  large  (Kp  = 
0->-6,  1-^5,  2-*-8,  5->-9)  . This  results  from  the  rigid  constraint  that 
all  points  to  be  lost  eure  lost.  Thus,  a point  which  just  escaped 
being  "peeled-off"  near  1500  LT  rotates  (at  a speed  less  than  the 
co-rotational  speed)  around  the  nightside  and  is  ultimately 
peeled-off  during  the  following  dayside  transit.  It  would  appear 
that  there  are  always  points  which  fall  into  this  category  (or 
linger  for  a very  long  time  near  the  stagnation  point  at  1800  LT) , 
and  thus  a diurnal  (=30  hr.)  time  scale  results  for  the  Tp(100%) 
condition. 

Waiting  for  the  "last  point  to  be  lost"  clearly  seems  to  be 
too  rigid  a condition.  Since  we  are  more  interested  in  the  plas- 
mapause  location  on  the  nightside  (where  its  relationship  to  the 
trough  may  best  exist) , and  since  the  nightside  sunward  convec- 
tion affects  the  L (18-<-06  LT)  values  more  directly,  we  view  the 
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50%  equilibriiim  values  of  Tp  to  be  more  representative  of  the 
time  lag  associated  with  a Kp/plasmapause  transition.  An  exam- 
ination of  the  TpC50%)  results  in  Table  2 reveals  that  Tp  becomes 
smaller  the  larger  the  change  in  Kp.  Thus,  when  AKp  = 1,  Tp  = 
11-16  hours,  for  AKp  = 2,  Tp  = 9-10  hours,  for  AKp  = 4,  Tp  = 6-8 
hours,  and  for  AKp  = 6,  Tp  = 4-6  hours.  This  seems  to  be  a con- 
sistent relationship  for  all  Kp  levels  and  may  be  described  in 
the  average  by  the  curve  given  in  Figure  5.  The  results  suggest 
that  if  statistical  studies  relating  the  nighttime  plasmapause 
to  Kp  are  to  be  used  as  predictors  of  actual  boundary  locations, 
an  appropriate  time  lag  should  be  used  dependent  upon  the  actual 
increase  in  Kp  encountered. 

2.3.  Simulation  of  Plasmapause  Disturbance  Scenarios'. 

The  second  aspect  of  our  investigations  using  Program  MAGCON 
was  an  attempt  to  simulate  actual  "disturbance  scenarios"  via  the 
model  cind  relate  the  changes  observed  to  actual  measurements 
CISIS-1  data  furnished  by  AFGL)  and  past  statistical  studies.  We 
selected  several  case  study  periods  for  simulation,  the  choices 
being  determined  by:  (1)  availability  of  AFGL  data,  (2)  the  de- 
sire to  have  well-defined  disturbance  periods  following  a period 
of  steady  activity  (i.e.,  Kp  = 1,2  or  3 for  nearly  24  hours),  and 
C3)  the  realization  that  only  the  2200-0600  LT  period  lends  itself 
to  our  modelling  capabilities. 

Before  proceeding  we  had  to  contend  with  the  fact  that  the 
equilibrium  plasmapause  contours  given  in  Figure  4 do  not  give 
"reasonable"  Lpp  values  in  the  pre-midnight  to  dawn  sector.  As 
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mentioned  earlier,  the  "tear-drop"  shape  gives  too  severe  a dawn- 

dusk  asymmetry.  The  50%  decrease  in  L from  1800  LT  to  00  LT  is 

PP 

not  reflected  in  actual  data,  and  thus  we  decided  to  re-normalize 

our  quadratic  vs.  Kp  relation  (as  derived  in  Mendillo  and 

Papagiannis  C1971)  for  the  1800  LT  sector)  to  actual  L (00  LT) 

PP 

AFGL/ISIS  observations.  An  examination  of  the  existing  data 

base  revealed  that  L (00  LT)  = 5 for  K =2,  and  thus  this 

PP  P 

normalization  point  was  chosen  for  the  E^^j^  vs.  Kp  relation,  rather 

than  the  Lpp  (18  LT)  = 5 for  Kp  = 3 used  previously.  The  new  set 

of  reduced  electric  field  values  and  consequent  L locations  are 

PP 

given  in  Table  3. 

To  illustrate  the  procedure,  consider  the  storm  period  of 
21-22  April,  1971.  Prior  to  1200  UT',  magnetic  activity  was 
steady  for  nearly  9 hours,  with  Kp  = 2.  From  1200  UT  on,  Kp  (a 
3-hour  index)  varied  as  follows:  3,  4,  5,  7,  5,  3,  2.  To  follow 
the  nightside  Lpp  variations  associated  with  this  scenario,  we 
first  used  our  simulation  program  to  generate  the  equilibrium 
C" tear-drop")  particle  distribution  appropriate  to  the  t=0  point 
of  the  scenario.  Thus,  the  initial  circular  distribution  of 
particles  was  allowed  to  evolve  to  the  case:  12:00  UT  on  21  April 
1971,  Kp  = 2 (Figure  6) . The  electric  field  model  was  then  used 
to  input  new  Ej^^  fields  in  concert  with  the  observed  Kp  variation. 
Computer  generated  plots  of  the  modified  particle  distributions 
were  made  at  hourly  intervals , and  the  Lpp  boundary  values  for 
the  LT  sector  in  question  were  scaled  from  the  plots  (see  Figure 
7 for  sample  outputs  2,  9,  12,  and  20  hours  into  the  simulation). 
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TABLE  3 

EQUILIBRIUM  PLASMAPAUSE  MODEL  NORMALIZED  TO  MIDNIGHT  SECTOR 


Kp 

V Csw) 

Ocm/sec) 

E 

conv 

Cmv/m)  Ckv/rel 

L 

PP 

Cdusk) 

L 

PP 

(noon/mi dn) 

Sp 

(dawn) 

0 

330.0 

0.07 

0.46 

14.09 

7.05 

5.84 

1 

397.5 

0.10 

0.66 

11.  70 

5.85 

4.85 

2 

465.0 

0.14 

0.91 

10.00 

5.00 

4.14 

3 

532.6 

0.19 

1.19 

8.  73 

4.  37 

3.62 

4 

600.1 

0.24 

1.51 

7.75 

3.  87 

3.21 

5 

667.6 

0.29 

1.  87 

6.97 

3.48 

2.87 

6 

735.1 

0.36 

2.27 

6.33 

3. 16 

2.62 

7 

802.6 

0.42 

2.70 

5.79 

2.90 

2. 40 

8 

870.2 

C.50 

3.18 

5.34 

2.67 

2.21 

9- 

937.7 

0.58 

3.69 

4.96 

2.48 

2.05 
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FIGURE  7. 


deformations  of  the  plasmasphere 
depicted  in  Figure  8 at  t=2,9,12  and  20  hours 
into  the  Kp  scenario  for  21-22  April  1971.  The  arrows 
indicate  the  Lpp  boundary  segments  as  defined  by  the 
particles  at  23:45  LT. 
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For  this  particular  scenario,  we  examined  the  2345  LT  values  in 
order  to  compare  the  model's  predictions  with  ISIS-1  in-situ 
probe  data  for  the  equatorward  edge  of  the  total  ion  trough  in 
the  topside  F-region  (Wildman  et  al.,  1976).  The  results  are 
given  in  Figure  8,  together  with  the  "statistical  predictions"  of 
the  L-value  of  the  plasmapause  using  the  relation  given  by 
Rycroft  and  Thomas  (1970)  . 

In  examining  the  results,  there  are  several  features  to  note: 
Cl)  The  trough's  equatorward  edge  (measured  high  in  the  F-region) 
is  probably  a better  indicator  of  the  plasmapause  location  during 
disturbed  times  than  during  quiet  periods.  Thus,  at  the  t=0  point 
of  the  scenario,  the  difference  between  the  data  and  the  model  is 
the  largest. 

(2)  As  the  storm  progresses  (in  this  case,  rather  smoothly  from 
Kp  = 3-*-4-^5-^8)  , the  steady  contraction  of  the  plasmapause  is  re- 
produced. The  statistical  results  of  Rycroft  and  Thomas  (1970) 
appear  as  discontinuous  jumps  every  3 hours  (the  Kp  interval); 
they  also  tend  to  agree  with  the  observed  behavior. 

(3)  Following  the  storm  maximvim,  i.e.,  after  00  LT  on  22  April 
1971,  the  Kp  transition  5->3->-2  causes  a gradual  recovery  of  the 
plasmapause  to  higher  L-values . There  is  a definite  "sluggish- 
ness" built  into  the  model  arising  from  the  previously  high  Kp 
values.  This  ability  to  take  previous  magnetic  history  into 
account  is  in  sharp  contrast  to  the  statistical  predictions  which 

show  instamtaneous  adjustments  to  higher  L values. 

PP 
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We  have  tested  this  procedure  for  other  Kp  scenarios  and 
find  the  "sluggishness"  of  the  model  to  be  a major  asset  for  both 
intensification  effects  (e.g.,  Kp  = l->-5)  and  recovery  patterns. 

Figure  9 gives  a second  example  for  which  frequent  satellite 
data  were  available.  The  period  considered  was  a large  and  pro- 
tracted storm,  giving  us  the  opportunity  to  test  the  modelling 
scheme  for  an  extended  period  of  relatively  high  geomagnetic 
activity.  The  model  predictions  are  again  remarkably  consistent 
with  the  satellite  data. 

In  summary,  our  numerical  simulation  studies  of  magneto- 
spheric  convection,  as  related  to  plasmapause  dynamics,  have 
yielded  some  quantitative  conclusions  about  plasmapause  transi- 
tion times  and  AKp  (summarised  in  Figure  5) , and  definite  evi- 
dence that  a relatively  simple  physical  model  for  a convection- 
dominated  plasmapause  can  be  used  to  simulate  nightside  dynamics 
effects  (as  evidenced  in  Figures  8 and  9) . 
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3.  FINITE  ELEMENT  SIMULATION  APPLIED  TO  THE  AURORAL  IONOSPHERE 

Two  of  the  major  effects  precipitating  magnetospheric  par- 
ticles have  upon  the  earth's  upper  atmosphere  are  excitation  of 
the  auroral  spectriom  and  the  production  of  ionization.  VVhile 
auroral  studies  have  grown  into  an  important  discipline  in  space 
research  Csee,  for  example,  the  treatise  by  Akasofu,  1977)  the 
latter  effect  has  attracted  only  limited  attention  [Rees,  1963;  Banks 
et  al.,1974;  Mantas  and  Walker,  1976]  despite  its  importance  in 
applications  to  radio  wave  propagation.  The  processes  control- 
ling the  redistribution  of  particle-generated  ionization  at  high 
latitudes  are  many  and  complex:  they  stem  from  chemical  reac- 
tions, vertical  diffusion,  heating,  E;- fie  Ids,  neutral  winds  and 
related  effects. 

The  purpose  of  the  investigation  reported  on  here  has  been 
to  develop,  as  a first  step  towards  a more  comprehensive  capa- 
bility, a one-dimensional  simulation  technique  to  describe  the 
vertical  distribution  of  electron  density  in  the  high  latitude 
ionosphere.  The  computational  procedure  adopted  is  Finite 
Element  Simulation  [Balko  and  Mendillo,  1977] . The  FES  technique 
relies  on  the  application  of  physical  and  chemical  laws,  usually 
in  a linearized  version,  to  the  system  under  investigation  in 
the  following  manner. 

The  system  is  divided  into  a number  of  convenient  cells 
Clabeled  i,j,k),  using  the  symmetry  of  the  situation  to  advantage 
so  as  to  minimize  the  number  of  geometrical  parameters. 
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Each  cell  is  characterized  by  three  types  of  parameters: 

a)  Geometrical  parameters  5(i) , such  as  the  position  of 
the  center  of  a cell/  its  volume,  cross  sectional  area, 
surface  area,  etc. 

b)  Physical  parameters  ^Ci) , such  as  diffusivity,  specific 
heat,  etc. 

c)  Transport  functions  rCi),  such  as  temperature,  concen- 
tration, pressure,  etc. 

T Ci)  is  an  intensive  parameter  and  the  associated  extensive 
parameters,  e.g.  heat,  number  of  particles,  is  denoted  Q(i). 

The  parameters  in  group  (c)  are  specified  initially,  t (i)  at 
t=0,  and  then  a time  interval  At  is  chosen  over  which  the  cal- 
culation is  to  take  place . The  choice  of  At  is  governed  by  the 
requirement  that  it  should  be  small  enough  to  prevent  "overshoot" 
in  the  calculations . 

The  heart  of  the  calculation  is  the  set  of  physical  laws 
which  describe  changes  in  the  system  parameters  T(i)  over  the 
interval  (t,  t + At) . A general  process  (m)  can  be  thought  of 
as  causing  a change  AQ  '^(i)  due  to  the  transfer  parameter  T”'(i) 
of  cell  i.  It  can  be  calculated  from 

AQ  ”'a)  = fid,  i + 1)  • A(i,  i + 1)  ’ $(i,  i + 1)  At  (1) 
where  fi,  A and  i are  functions  of  the  appropriate  geometrical, 
physical  and  transfer  parameters  of  the  system. 

AQ  ”*Ci)  is  a change  in  the  extensive  parameter  corresponding 
to  T^^d)  and  gives  the  amount  of  "material"  transferred  out  of 
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cell  Cil  by  the  process  m.  In  the  case  of  molecular  diffusion, 

AQ  is  the  number  of  molecules  transported  in  time  At,  due  to 
concentration  differences  in  the  cells.  The  changes  in  the 
transport  functions  can  be  obtained  from 

AT"‘Ci)  = CAQ"‘(i-l)  - AQ^^Ci)  )/R"*(i)  (2) 

This  leads  to  a new  t Cil  at  time  t + At  which  is  obtained  from 
all  the  At™ Cil  by 

ra=l 

where  we  Rave  assrnned  OH  independent  processes  affecting  the 
parameter  r CH • 

The  preceding  three  steps  are  repeated  until  is  reached 

The  result  of  the  calculation  is  a profile  of  the  x parameter  as 
a fvmction  of  position  (cell  center)  for  different  times.  Note 
that  this  is  a dynamic  simulation  calculation  so  that  in  order 
to  obtain  an  equilibrium  situation  the  temporal  development  has 
to  be  simulated  until  an  asymptotic  behavior  is  reached. 

One  Dimensional  Diffusion 

Figure  10  represents  the  simple  geometry  we  have  adopted 
for  the  problem  of  vertical  diffusion  of  charged  peirticles  in 
the  earth's  upper  atmosphere  C150  ^ h _<  1500  km).  For  simplic- 
ity  we  assume  1)  the  earth's  magnetic  field  B = BZ  and  2)  the 
horizontal  cross-sectional  curea  of  the  cells  A(Z)  = A = constant 
Typically,  A =»  20  km  x 20  km  and  the  total  number  of  cells  = 40. 
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Figure  10,  Cell  Geometry  for  One-Dimensional  Vertical 
Diffusion. 
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The  number  of  molecules  transferred  by  diffusion  from  cell 
Ci-1)  to  cell  Ci)  during  the  time  interval  At  is  given  by 

ANCi-1,  i)  = [C(i)  - C(i-l)  ] • At  (4) 

AZCi-1/  i) 

where  for  cell  i the  quantity  ANCi-l»  i)  is  positive  if  the  par- 
ticle flux  is  into  the  cell.  Note  that  the  diffusivity  D repre- 
sents effective  diffusivity  defined  as 

i.  T 

^ (5) 

° - H/Ol 

where  Ai  is  the  length  of  the  i^  cell.  C(i)  is  the  concentration 
of  the  diffusing  species  in  cell  i,  AZ(i-l,  i)  = Z (i)  - Z(i-l), 
separation  of  centers  of  cells  i and  i-1. 

Analytically,  this  expression  can  be  written  as  a flux  from 
one  cell  into  the  other,  dependent  on  the  gradient  between  the 
two  cells.  Defining  a flux  into  a cell  as  positive, 

F = -D  dC 

az 

where  D is  again  the  diffusion  coefficient  and  F is  the  number 

of  particles  crossing  a unit  area  per  unit  time  between  the 

interface  of  the  two  cells,  F = ^ dN  . The  concentration  change 

A at 

in  the  ith-cell  after  diffusive  transfer  between  it  and  the  two 
adjacent  cells  i-1  and  i+1  is 


AC°  = IAN°Ci-l,  i)  - AN°Ci,  i+D]  • — (( 

^i 

where  VCi)  is  the  volume  of  the  ith  cell.  This  is  a special 
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case  of  equation  (2)  where  the  parameter  is  simply  the 

volume  of  the  ith  cell.  The  new  concentrations  at  time  t + At 
are 


C^(t  + At)  = C^(t)  + AC?  (7) 

This  operation  is  performed  for  all  the  N cells  in  the  system 

and  is  repeated  n times  for  s\ibsequent  time  intervals  until  the 

maximvim  required  time  = nAt  is  reached.  The  calculation 

step  CAt)  need  not  remain  constant  over  the  entire  time  span 

could  be  changed  as  the  calculation  evolves.  The  re- 
mcoc  ^ 

suit  of  the  calculation  is  a "concentration  profile"  of  the 
system  as  a function  of  time. 

We  can  now  show  the  correspondence  between  the  above  approach 
and  the  diffusion  equation.  Consider  the  simple  case  of  a 
homogeneous  medium.  The  concentration  change  over  the  time 
interval  At  as  calculated  from  equations  (4)  and  (6)  may  be 
written  as 


AC^Ct,  At) 
At 


= ■=  • s|-'=i-i  * - 2Cil  • k: 


(8) 


where  = A * AZ^  is  the  volume  of  cell  i and  all  cells  are 
assumed  to  have  identical  cross-sectional  areas  and  linear 
dimensions.  Changing  notations  to  C^  C(Z),  -*■  C(Z-AZ), 

-*■  C (Z  + AZ)  and  taking  the  limit  as  At  -►  0 and  AZ  -*•  0 
independently,  we  get 


D 


2 


3 C 


l£ 

3 t 


(9) 
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by  the  definition  of  the  derivative.  We  have  thus  derived 
the  diffusion  equation  from  our  two  equations  of  transport  in  a 
hcsnogeneous  medium. 

The  next  step  in  building  up  our  simulation  of  the  trans- 

F 

port  phemomena  is  to  introduce  the  force  field  drift  term  (ACj^) 
which,  in  our  case  of  atmospheric  diffusion,  arises  from  the 
Earth's  gravitational  attraction.  Given  an  effective  drift 

I 

velocity  between  cells  (due,  in  this  case,  to  the  gravita- 
tional acceleration)  the  number  of  molecules  transferred  from 
cell  Ci)  to  cell  (i+1)  in  the  time  At  is  given  by 

P . 

AN^Ci,i  + 1)  + G. — ^C.  At 

AZi 

where  G^  is  calculated  the  same  way  as  (see  equation  5) . The 
change  in  the  concentration  of  cell  i due  to  this  process  is 

AC^  = [AN^(i-l,  i)  - AN^(i,  i +1)]  • ^ (10) 

'^i 


Analytically,  this  change  is  usually  referred  to  as  the  ' advection 
term'  in  the  continuity  equation  for  concentration  (C) . It  gives 
the  time  rate  of  change  in  C due  to  drifting  spatial  gradients 
in  C,  i.e. 


3 t 


V • VC 


3Z 


Following  the  procedure  used  with  the  diffusive  term  above,  we 
can  show  that  equations  (1) , (2) , (5) , (10) , in  the  appropriate 
limit  can  be  used  to  derive  Smoluchowski' s equation  for  diffu- 
sion in  a force  field 
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iC  ^ D 31^2+  G 
3t  3Z  3Z 


(11) 


for  a homogeneous  medium. 

Finally,  for  reactive  substances  we  introduce  the  change 
in  concentration  due  to  chemical  reactions.  Suppose  the  diffus- 
ing substance  reacts  with  a background  constituent  whose  concen- 
tration  in  cell  i is  given  by  . The  change  in  concentration 
of  due  to  reactions  of  the  diffusing  material  during  the  time 
interval  fit  is  given  by 

fiC?  = -kC.  • fit  (12) 

1 11 

where  k is  the  chemical  reaction  rate  in  units  of  volume/time. 

In  summary,  if  Ci*  represents  the  concentration  of  substance 
« C * = 1, . . .n ^ ) in  cell  i,  then  the  total  change  in  concentra- 
tion is  obtained  from 

fiCi  (t,  fit)  = fiqf,i  + fiqf,i  + Aqf,i  (13) 

Computations  and  Results 

Electron  density  enhamcements  in  the  F-region  beneath  the 
day  side  magnetospheric  cusp  represents  a relatively  simple  case 
of  particle-generated  ionization  in  the  earth's  upper  atmosphere 
{Chacko  and  Mendillo,  1977] . We  have  applied  the  FES  technique 
sketched  in  the  previous  section  to  this  problem.  The  model 
incorporates  vertical  diffusion  of  o''’  and  H*^,  gravitational 
drift  and  predominant  chemical  reactions. 
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Table  4 shows  values  of  the  diffusion  coefficients  we  have 
adopted  iBanks  and  Kockarts,  1973] . Table  5 lists  the  chemical 
reactions  and  their  rate  coefficients  included  in  the  computa- 
tions . 

Figure  11  shows  the  initial  and  distributions  which 
form  part  of  the  input  into  the  program.  The  lower  parts  of 
both  profiles  as  well  as  the  upper  part  of  the  profile  re- 
present reasonable  but  arbitrary  values  while  the  upper  F-region 
values  for  [0^]  are  the  measured  electron  densities  at  75.5° 
Corrected  Geomagnetic  Latitude  during  quiet  periods  in  December 
1971  [Chacko  and  MendillO/  1977] . These  initial  profiles  are 
considered  adequate  for  the  purpose  of  demonstrating  the  feasi- 
bility of  the  FES  scheme  in  the  area  of  plasma  diffusion  at  high 
latitudes . 

Figxire  12  shows  the  evolution  of  0^  and  h'*’  profiles  at  the 
end  of  15  minutes  as  a result  of  diffusion  and  the  chemical 
reactions  listed  above  as  well  as  of  gravitational  drift. 

Figure  13  represents  the  vertical  0^  and  distributions  at 
the  end  of  20  minutes  of  real  time  when,  in  addition  to  the 
cibove  processes,  in-situ  production  by  soft  electrons  [Knudsen 
et  al.,  1977]  was  allowed  to  operate  for  the  last  5 minutes. 

The  fact  that  the  final  profiles  are  well-behaved  and  quantita- 
tively realistic  shows  that  the  FES  scheme  we  have  developed 
CeUi  now  be  utilized  to  obtain  quantitative  results  in  the  field 
of  plasma  diffusion  at  high  latitudes. 
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TABLE  4 


DIFFUSION  COEFFICIENTS 


Species 

D 

2 

Cc3U  /sec) 

0 
+ 

1 

o 

3.1  X 

10^^ 

T 

[OJ 

(1  + t^/t^)  ^ 

0 
+ 

1 

o 

7.2  X 

10 

T . 

O'^-Na 

8.5  X 

10^5 

T. 

TnIT 

0'''-He 

5.6  X 

10 

TSiT 

0^-H 

4.0  X 

10 

ThT 

0 
+ 

1 

+ 

9.4  X 

10^ 

T 5/2 

1 

[H^] 

+ 

1 

o 

3.7  X 

10 16 

I’i 

+ 

1 

o 

2.6  X 

10l6 

H'''-N2 

2.5  X 

I0l6 

TnIT 

H'''-He 

7.8  X 

10l6 

TijiT 

+ 

1 

o 

+ 

6.9  X 

10^ 

t5/2 

i 

[o'"] 

h‘'‘-h 

8.0  X 

loi^ 

’■i'’ 

IHJ 

(1  + T^/T^)^ 
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TABLE  5 


CHEMICAL  REACTIONS  AND  THEIR  RATES 


o'*’  + N2  n + no'*’ 

0 + O2  -►  0 + O2 

no'*’  + e“  -►  N + 0 

Oz^  + e -*-0  + 0 

o'*”  + H h'*’  + 0 

e'*’  + 0 o'*'  + H 


1.2  X 10"^^(300/T^) 

2.0  X 10"^^(300/Tj^)  ^ 

4.1  X 10"^(300/T^) 

2.2  X 10"'^(300/Tg)°*®^ 

2.5  X 10"^^  T ^ 
n 

2.3  X 10“^^  T^^ 
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(m^) 


Figure  11.  The  initial  (t  = 0)  ionospheric  density  profiles.  The  electron 

densities  represent  measured  values  at  75.5*  Corrected  Geomagnetic 
Latitude  on  the  northern  day  side  in  midwinter  during  quiet  magnetic 
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